Chapter 8 Captive-bred vs wild

load("data/beta_filtered_capwild.Rdata")
wild_pre_samples <- sample_metadata %>% 
  filter(diet!="Post_grass") %>%
  dplyr::select(sample) %>% pull()

genome_counts<- genome_counts %>% 
  column_to_rownames("genome") %>% 
  select(all_of(wild_pre_samples))%>% 
  rownames_to_column("genome")

sample_metadata <- sample_metadata %>% 
  filter(diet!="Post_grass")
beta_q0n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

genome_counts <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)

beta_q1p <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt %>%
#  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

8.1 Permanova analysis

set.seed(1234)

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.004224 0.0042240 7.8055    999  0.003 **
Residuals 22 0.011905 0.0005412                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.008
Wild       0.010583      
adonis2(beta_q0n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_y2dqa26z39giqc2kfisq
term df SumOfSqs R2 statistic p.value
diet 1 0.01526961 0.2489907 7.29391 0.001
Residual 22 0.04605641 0.7510093 NA NA
Total 23 0.06132602 1.0000000 NA NA

Neutral

betadisper(beta_q1n$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.006925 0.0069254 0.5404    999  0.491
Residuals 22 0.281957 0.0128162                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.476
Wild        0.47005      
adonis2(beta_q1n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_mlxrdgkl0941ux9o1nm6
term df SumOfSqs R2 statistic p.value
diet 1 1.538238 0.2539205 7.487472 0.001
Residual 22 4.519715 0.7460795 NA NA
Total 23 6.057953 1.0000000 NA NA

Phylogenetic

betadisper(beta_q1p$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.007671 0.0076708 1.8217    999  0.179
Residuals 22 0.092637 0.0042108                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.187
Wild        0.19084      
adonis2(beta_q1p$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_3f450bngrbc1njb7kdqe
term df SumOfSqs R2 statistic p.value
diet 1 0.3720264 0.3082818 9.804858 0.001
Residual 22 0.8347476 0.6917182 NA NA
Total 23 1.2067741 1.0000000 NA NA

Functional

betadisper(beta_q1f$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.003070 0.0030698 1.7361    999  0.171
Residuals 22 0.038901 0.0017682                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.176
Wild        0.20119      
adonis2(beta_q1f$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_glmf4taaiuawxafp2s1w
term df SumOfSqs R2 statistic p.value
diet 1 0.01079711 0.187686 5.083124 0.103
Residual 22 0.04673041 0.812314 NA NA
Total 23 0.05752752 1.000000 NA NA

8.2 Beta diversity plots

8.2.1 Richness

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet", "Wild" = "Wild"))+
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(face = "bold", size = 12),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 12),
      legend.title = element_text(size = 14),
      legend.position = "right", legend.box = "vertical"
    ) +labs(color='Origin')

8.2.2 Neutral

8.2.3 Phylogenetic

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

8.2.4 Functional

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
  geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme_classic() +
  theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')